#include "qpcrModel.h"

using namespace QPCRModel;

#pragma region CONSTRUCTOR{

Model::Model()
{
    // 初始化参数值
    _cycleNum = -1;
    _RNoise = 0;
    _CSPE = -1;
    _CFDM = -1;
    _CSDM = -1;
    _E = 0;
    _slope = 0;
    _Rmax = 0;
    _Rmin = 0;
    _fittingConf = -1;
    _FINISHED = false;
    _S_SHAPED = false;

    _functype = 1;
    _expLen = 6;
    _baseLen = 3;

    _pSigmoid = new funcSigmoid;
    _pRichards = new funcRichards;
    _pLogistic = new funcLogistic;
    _pExponent = new funcExponent;
    func = _pRichards;

}

Model::~Model(){

    if(_pSigmoid != nullptr){
        delete _pSigmoid;
    }
    if(_pRichards != nullptr){
        delete _pRichards;
    }
    if(_pLogistic != nullptr){
        delete _pLogistic;
    }
    if(_pExponent != nullptr){
        delete _pExponent;
    }
}

#pragma endregion CONSTRUCTOR}


#pragma region INTERFACE{


void Model::setRawData(QVector<QPointF> pointList){

    std::sort(pointList.begin(), pointList.end(), [](QPointF a, QPointF b){return a.x() < b.x();});

    _rawPointList.clear();
    int c = 0;
    for(auto it = pointList.begin(); it != pointList.end(); ){
        if(c < it->x()){
            _rawPointList.append(QPointF(c, it->y()));
            c ++;
        }else if(c > it->x()){
            it ++;
        }
        else{
            _rawPointList.append(QPointF(c, it->y()));
            c ++;
            it ++;
        }
    }

    _cycleNum = _rawPointList.size()-1;
    if(_cycleNum < 10){
        //qDebug("ERROR - Sigmoid::setRawData: 数据输入过少");
        _cycleNum = -1;
        return;
    }

}

void Model::setParam(ModelParam type, double value){
    if(type == FUNC_TYPE){
        if(value == 0){
            func = _pSigmoid;
            _functype = value;
        }else if(value == 1){
            func = _pRichards;
            _functype = value;
        }else if(value == 2){
            func = _pLogistic;
            _functype = value;
        }
//        else if(value == 3){      // 指数函数拟合暂不开放
//            func = _pExponent;
//            _functype = value;
//        }
    }
    else if(type == BASE_LEN){
        value = std::max(value, 1.0);
        value = std::min(value, 10.0);
        _baseLen = value;
    }


    else if(type == BASE_START){
        value = std::max(value, 0.0);
    }else if(type == BASE_END){
        value = std::max(value, 0.0);
    }

}

void Model::startCalculation(){

    if(_cycleNum < 10){
        return;
    }

    _FINISHED = false;

    _S_SHAPED = calculate() >= 0 ? true : false;
    updatePointList();

    _FINISHED = true;
}

QVector<QPointF> Model::getRawData(){
    return _rawPointList;
}

QVector<QPointF> Model::getFittedData(){
    //return filter(filter(filter(filter(_rawPointList))));
    return _fittedPointList;
}

QVector<QPointF> Model::getFittedData_delta(){
    return _fittedPointList_delta;
}


QVector<QPointF> Model::getFittedData_log(){
    return _fittedPointList_log;
}

QVector<QPointF> Model::getFittedData_log_delta(){
    return _fittedPointList_log_delta;
}

QVector<QPointF> Model::getFittedData_norm(){
    return _fittedPointList_norm;
}

QVector<QPointF> Model::getFittedDate_diff(double range){
    QVector<QPointF> res = _fittedPointList_delta;
    for(int i = 0; i < res.size(); i ++){
        if(range == 0){
            res[i].setY(0);
        }else{
            res[i].setY( res[i].y() / range );
        }

    }
    return res;
}

QVector<QPointF> Model::getBaselineData(){
    return _baselinePointList;
}

double Model::getSPE(){
    if(_fittingConf <= 0) return 0;
    return _CSPE;
}

double Model::getSDM(){
    if(_fittingConf <= 0) return 0;
    return _CSDM;
}

double Model::getFDM(){
    if(_fittingConf <= 0) return 0;
    return _CFDM;
}

double Model::getCt(){
    if(_fittingConf <= 0) return 0;
    return func->g(_threshold);
}


double Model::getCt(double threshold, YAxisType type){

    if(_fittingConf <= 0) return 0;

    switch(type){

    case RN:
        break;
    case DELTA_RN:
        threshold += _bgValue;
        break;
    case LOG_RN:
        threshold = pow(10, threshold);
        break;
    case LOG_DELTA_RN:
        threshold = pow(10, threshold) + _bgValue;
        break;
    case NORM:
        threshold = threshold * (func->f(_cycleNum)-func->f(1)) + func->f(1);
        break;
    }

    return func->g(threshold);
}



double Model::getCt_relative(int type){
    if(_fittingConf <= 0) return 0;
    /* TYPE 0: 0.5 * (F(CSDM) + F0)
     * TYPE 1: CSDM
     * TYPE 2: Cy0
     */

    if(type > 2 || type < 0){
        type = 0;
    }

    if(type == 0){
        double threshold = 0.5 * func->f(_CSDM) + 0.5 * _bgValue;   // 相对阈值
        return func->g(threshold);
    }else if(type == 1){
        return _CSPE;
    }else if(type == 2){    // Cy0
        double c = func->x0, b = func->b, d = func->d, Fb = func->y0, Fmax = func->a+func->y0;
        double lk, lb;
        func->getLinearStageFunc(&lk, &lb);
        double cy0 = c + b*log(d) - b * ((d+1)/d);
        cy0 = (func->y0 - lb) / lk;
        return cy0;
    }
    return 0;
}

double Model::getBgValue(){
    if(_fittingConf <= 0) return _Rmin;
    return _bgValue;
}

double Model::getThreshold(YAxisType type){
    if(_fittingConf <= 0) return 0;

    if(type == RN){
        return _threshold;
    }else if(type == DELTA_RN){
        return _threshold - _bgValue;
    }else if(type == LOG_DELTA_RN){
        return log10(_threshold - _bgValue);
    }else if(type == LOG_RN){
        return log10(_threshold);
    }else if(type == NORM){
        return (func->f(_cycleNum)-_threshold) / (func->f(_cycleNum)-func->f(1));
    }

    return _threshold;
}


double Model::getThreshold(int cycle, YAxisType type){
    if(_fittingConf <= 0) return 0;

    double height = func->f(cycle) - _bgValue;


    if(_fittingConf <= 0) return 0;

    if(type == RN){
        return _threshold + height;
    }else if(type == DELTA_RN){
        return _threshold - _bgValue + height;
    }else if(type == LOG_DELTA_RN){
        return log10(_threshold - _bgValue + height);
    }else if(type == LOG_RN){
        return log10(_threshold + height);
    }else if(type == NORM){
        return (func->f(_cycleNum)-_threshold) / (func->f(_cycleNum)-func->f(1));
    }

    return _threshold;
}


double Model::getEfficiency(){
    if(_fittingConf <= 0) return 0;
    return _E;
}

double Model::getSlope(){
    if(_fittingConf <= 0) return 0;
    return _slope;
}

QVector<double> Model::getFuncParam(){
    QVector<double> res;
    res.append(_functype);
    res.append(_fittingConf);
    res.append(func->a);
    res.append(func->b);
    res.append(func->x0);
    res.append(func->y0);
    res.append(func->d);
    return res;
}

void Model::printParam(){
    qDebug("PARAMETER:");
    qDebug("  SPE|SDM|FDM:  %.2f %.2f %.2f", _CSPE, _CSDM, _CFDM);
    qDebug("  bgValue:      %.2f", _bgValue);
    qDebug("  threshold:    %.2f", _threshold);
    qDebug("  ct_absolute:  %.2f", getCt());
    qDebug("  ct_relative:  %.2f | %.2f | %.2f", getCt_relative(0), getCt_relative(1), getCt_relative(2));
    qDebug("  E|Slope:      %.2f | %.2f", _E, _slope);
    qDebug("  FuncParam:    a=%.2f b=%.2f x0=%.2f y0=%.2f d=%.2f", func->a, func->b, func->x0, func->y0, func->d);
    qDebug("  Confidence:   %.2f", _fittingConf);

    // 原始数据
//    for(int i = 0; i <=_cycleNum; i ++){
//        qDebug("%.3f", _rawPointList[i].y());
//    }

//    QVector<double> funcParam = getFuncParam();
//    qDebug("staticCt: %f", Model::getCt(_threshold,funcParam, RN));

//    qDebug("  Thresholds:   %.2f %.2f %.2f %.2f", getThreshold(RN), getThreshold(DELTA_RN), getThreshold(LOG_RN), getThreshold(LOG_DELTA_RN));
//    qDebug("  setThreshold: A-%.2f B-%.2f C-%.2f D-%.2f", getCt(12000, RN), getCt(1000, DELTA_RN), getCt(4.1, LOG_RN), getCt(2.82, LOG_DELTA_RN));

    qDebug("FINDEXP:");
    for(int i = 0; i <= _cycleNum-5; i ++){
        qDebug("  cycle%d K = %.3f", i, _slopeList[i]);
    }
    return;
}

bool Model::finished(){
    return true;
}

double Model::getFittingConf(){
    return _fittingConf;
}

double Model::getTopValue(int type){

    if(_fittingConf <= 0) return _Rmax;

    if(type == 0){
        return func->f(_cycleNum);
    }else if(type == 1){
        return func->a + func->y0;
    }
    return func->f(_cycleNum);
}

#pragma endregion}
